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In materials science the phase field crystal approach has become popular to model crystallization 
processes. Phase field crystal models are in essence Landau-Ginzburg-type models, which should be 
derivable from the underlying microscopic description of the system in question. We present a study 
on classical density functional theory in three stages of approximation leading to a specific phase 
field crystal model, and we discuss the limits of applicability of the models that result from these 
approximations. As a test system we have chosen the three-dimensional suspension of monodisperse 
hard spheres. 

The levels of density functional theory that we discuss are fundamental measure theory, a second- 
order Taylor expansion thereof, and a minimal phase-field crystal model. We have computed coex- 
istence densities, vacancy concentrations in the crystalline phase, interfacial tensions and interfacial 
order parameter profiles, and we compare these quantities to simulation results. We also suggest a 
procedure to fit the free parameters of the phase field crystal model. 

In brief, we conclude that fundamental measure theory is very accurate and can serve as a bench- 
mark for the other theories. Taylor expansion strongly affects free energies, surface tensions and 
vacancy concentrations. Furthermore it is phenomenologically misleading to interpret the phase 
field crystal model as stemming directly from Taylor-expanded density functional theory. 

PACS numbers: 82.70Dd,61.50Ah,71.15Mb 
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I. INTRODUCTION 



In materials science, the modelling of dynamic processes involving the growth of solid phases in melts or in the 
environment of another solid phase has been advanced using phase field models in the past years [3, HI- Here, the 
phase field ip is associated with an order parameter field that distinguishes between a solid and a liquid phase, and it 
is usually coupled with a density or concentration field g for that phase. The dynamics of g is conserved, following 
the equation 

§ = v-(r e j e ) (i) 

Here, the density current j e is the gradient of a chemical potential function which is assumed to be derivable from a 
free energy functional J- . T e is a mobility. In contrast to this, there is no conservation law for the order parameter 
and thus one can generically assume a nonconserved dynamic evolution of tp of the form 



dip - 5F[g, p] 

at v 8p(v) ■ [ ' 



In order to briefly explain the approach, we consider a one-component system able to form one fluid and one solid 
phase. The simplest free energy functional which gives us phase coexistence associated with smoothly varying profiles 
for g and ip across the phase boundary follows from a gradient expansion in the specific truncation 

F[0, <f] = j^j; J d 3 r (c 2 (Vp) 2 + f( Q , ip)) . (4) 

Here it is assumed that any inhomogeneity costs free energy through the gradient term in the order parameter field. 
Since there is no corresponding gradient term in the density, the free energy penalty corresponding to a change 
in the density field must be small, and this appears to be only possible if g is suitably coarse-grained from the 
microscopic density field. Consequently, the variations in the microscopic density field relevant for the free energy 
must be contained in the phase field tp which in turn should be derivable from the microscopic density through another 
coarse-graining procedure. We will return to that point below. The potential function /(p, ip) contains a double-well 
type expression with minima at tp = — 1 (fluid) and tp = 1 (solid), modified such that minimization with respect to g 
gives the input fluid and crystal coexistence densities gn{T) and g CI (T) which in general depend on the temperature 
T. A possible form is Q 

f(g, tp) = gi (tp) + i(l + gi{ip))UM + \{l - 92(<p))Mq) , (5) 

with 9l (tp) = {ip + lf{ip-l)\ (6) 
and g 2 (tp = ±1) = ±1 , g' 2 (ip = ±1) = . (7) 

Thus, the phase field approach is nothing but a slighly rewritten Landau-Ginzburg model for the fluid-solid phase 
transition. The formulation accomodates an empirical free energy density for the fluid phase (fn(g)) and the crystal 
phase (/ C r(f?)) yielding the required input coexistence densities. In the form of Eq. Q, the free energy contains the 
parameters F related to is the free energy scale in units of the thermal energy ksT (should be adjusted to the bulk 
free energy difference of solid and liquid) and the constant c 2 which can be fixed through the value of the liquid-solid 
surface tension. 

In such a way, nucleation and growth in simple systems can be addressed without resolving the details of the free 
energy for inhomogeneous systems. For the widely used hard sphere reference system (which will be examined more 
in detail in this work), this has been done in e.g. Ref. [||. 

Such a minimal Landau-Ginzburg description can be extended to more complex systems. For each new material 
component of the sytem, one needs to introduce a corresponding density field, and new phase fields for the fluid 
and solid phases, even when the solid phases just differ by their crystalline orientation. Thus the number of free 
energy and surface tension parameters quickly grows when the complexity of the system is increased. Even for the 
one-component system, empirical information on the anisotropy of the surface tension for different crystal faces in 
equilibrium with the fluid needs to be taken into account to set the parameters. Hence an important question is 
whether the phase field itself can be consistently treated in terms of the density field which stems from a microscopic 
foundation. 

In this paper, we investigate three formulations (or approximations) to classical density functional theory which deal 
with the microscopic particle density field and thus, in principle, constitute the underlying theoretical framework from 
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which a consistent phase field crystal description should arise. In our explicit calculations, we examine the coexistence 
properties, surface tensions and interface density modes in the hard sphere system. Understanding surface tensions 
for different interface orientations and the associated surface structure are important prerequisites for further studies. 
There are three main reasons to choose the hard sphere system for method comparison: (1) availability of very precise 
density functionals (fundamental measure theory, our first formulation). (2) empirical evidence that the crystal- 
liquid surface tensions of fcc(hcp) forming metals are largely of entropic origin and thus the packing of impenetrable 
cores plays an important role for the surface structure of these metals Q and (3) the athermal nature of the hard 
sphere system which reduces the parameters for describing the coexistence to just the pair of coexistence densities 
for the liquid and solid phase. - Our second formulation is Taylor-expanded density functional theory which neglects 
the density fluctuations with respect to a reference density beyond second order in the free energy formulation. It 
constitutes already a drastic approximation in density functional theory, nevertheless it is occasionally depicted in 
the literature as "the" density functional theory from which the third formulation investigated here, the phase field 
crystal model of the simplest type, can be derived. In short, the phase field crystal model can be viewed as a local 
expansion in density fluctuations and in their gradients of Taylor-expanded density functional theory. The model 
(see below for its description and references) has lately come to some prominence in the materials science community 
mainly for the reason alluded to a microscopic foundation of phase field descriptions, (see Ref. [BJ for a systematic 
attempt in that direction). 

The paper is structured as follows. Sec. |TT] introduces briefly the foundations of classical density functional theory 
and describes the formal approximation steps leading to Taylor-expanded density functionals and phase field crystal 
models. In Sec. IIII1 the explicit functionals are given and they are applied to the hard sphere system. In calculating 
surface tensions and the interface structure, no further approximations are made in order to avoid uncertainties in 
interpreting the results. We concisely discuss the problem of parameter fixing for the phase field crystal model. Sec. lIVI 
contains our summary and conclusions. 



II. DENSITY FUNCTIONAL THEORY AND PHASE FIELD CRYSTAL MODELS 



As discussed before, the phase and the density field in the phase field formulation should be both obtainable through 
a suitable coarse graining of the microscopic density field. Thus one may be tempted to forego the artificial distinction 
between (coarse-grained) phase and (coarse-grained) density field entirely in favor of the microscopic density p(r). 
Indeed, in equilibrium the basic theorems of density functional theory assure us that there is a unique free energy 
functional of the one-particle density field p(r), 

T[p}=T id [p}+T cx [ P }, (8) 
with pF id [p] = J d 3 rp(r) (ln(p(r)A 3 ) - l) (9) 

which can be split into the exactly known ideal gas part J 71 ' 1 (A is the de-Broglie wavelength, f3 = \/{k-QT) is the 
inverse temperature) and a generally unknown excess part J 76 *. The equilibrium density p oq in the presence of an 
external (one-particle) potential V cxt (r) is then given by 

6T[p] 



5p(r) 



(10) 



P = Poq 



where p is the imposed chemical potential (e.g. by requiring a certain bulk density far away from the region where 
the external potential acts). For diffusive dynamics, the time evolution of this microscopic one-particle density field 
obeys the type of dynamics as in Eqs. (fT|) and ©: 



p(r,t)V 



+ V cxt (r,t) 



(11) 



_Sp(r,t) 

To show this, one needs the possibly severe approximation that the time-dependent density-density correlation 
function can be approximated by the corresponding equilibrium object which in turn is obtainable from the equilibrium 
density functional F[p] Note that the density field p(r, t) is an ensemble-averaged quantity with no coarse-graining 
in space and time and there is no noise term in Eq. pip . 



A. Functional Taylor expansion 

Since the excess free energy functional J-" cx is unknown in general, many practical applications of DFT have started 
from an expansion of J rcx around a background reference density profile po( r ) which, in general, can depend on the 
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position: 

f3F c * = /3F ex [p ] - J (Pre™ ( r; p ) Ap(r) - ± / dVdVc® (r, r'; p ) Ap(r) Ap(r>) + ... (12) 

Here, Fq*[pq] is the excess free energy pertaining to the background profile, Ap(r) = p(r) — po(r) and c^ and c*- 2 -* 
are the first two members in the hierarchy of direct correlation functions , defined by 

(13) 

P=Po(r) 

In most practical applications, p = const, is taken to be a reference bulk density in which case — = /3/i cx = 
/3/i — log(p A 3 ) and c^ 2 \r — r'; p ) depends only on one position. To evaluate the functional in Eq. (fT2"j) . the correlation 
function c*- 2 -* has to be determined as an external input, provided e.g. by integral equation theory or by simple 
approximations of RPA type Q . 

It is perhaps somewhat surprising that the functional in Eq. (|12| (with po being a bulk density) is capable of 
describing fluid-solid coexistence. This has been shown first in Ref. [8[ for the case of hard spheres (for an fee crystal 
structure) with taken to be the analytically known solution of the Percus-Yevick closure to the integral equations. 
After all, the direct correlation function in the solid phase should be very distinct from the one in the liquid phase, as 
can be inferred from the definition in Eq. (|13[) . Consequently the expansion should only hold for modest departures 
from the reference bulk density which is not the case when comparing the density distribution in the crystal and 
the liquid, owing to the occurence of sharply peaked crystal density profiles. However, the fee crystal density profile 
probes the Fourier transform c^ 2 \k\ p) at discrete points in fc -space (the reciprocal lattice vectors (RLV)) which 
include the points where the structure factor has its maxima. Furthermore, the k-vectors of the RLV are distributed 
fairly isotropically (see also a more detailed discussion on that in Ref. Q). 

With suitable input for cS 2 \ the Taylor expanded functional in Eq. (fl"2|) is also capable of describing the fluid-bec 
transition (relevant for the description of e.g. iron). See Ref. [ic| for a recent study. However, the numerical results for 
free energies and also for crystal-liquid surface tensions obtained in such studies do not compare well with simulation 
results, e.g., the surface tensions from Ref. [To| deviate by a factor of 2. Thus, the Taylor-expanded functional appears 
to be merely a suitable qualitative tool to explore basic features of dense liquids in the vicinity of the solid or glass 
transition (see e.g. [TTI. Il2|). 
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B. The phase field crystal model 

The Taylor expanded functional in Eq. (|12p is nonlocal in the densities. Through an additional approximation 
(gradient expansion) it can be cast into a local form. We consider again a constant reference density po and the 
following power expansion of the Fourier transform of the direct correlation function: 

(k; p) = -c + c 2 k 2 - c 4 fc 4 . . . (14) 

Using this, the Taylor-expanded functional becomes 

/?Ji c x c = PF^{po) + f3p cx J d\A P {v) d 3 rAp(r) (c + c 2 V 2 + c 4 V 4 . . . ) Ap(r) + . . . 

(15) 

We observe that the excess free energy density contains local terms up to order 2 in Ap and up to order 4 in V(Ap). 
The total free energy contains additionally the ideal gas term, F ld [p] from Eq. ([9]). One may expand also this term 
in Ap in order to obtain a consistently power-expanded free energy density. It has been customary in the literature 
to introduce the dimcnsionlcss density difference as an order parameter: 

</>(r) = . (16) 

Po 

In terms of the power-expanded total free energy up to order 4 in and V</> reads 

I3T w f3F (po) + / d 6 r \p /3p,(/)(r) + + 

f / d 3 r</)(r) (c + c 2 V 2 + c 4 V 4 ) 0(r) . (17) 
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The coefficients a; d = p , = Po/2, and g 1( \ = p /3 stem from the expansion of the ideal gas part of the free 
energy. The reference free energy F contains the ideal gas part by F (p ) = F^ x (p ) + F ld (p ). The model defined 
in Eq. ([T7| looks like a straightforward extension of standard square-gradient Ginzburg-Landau models. It has been 



formulated in Ref. (l3|, however, earlier work has established the usefulness of such a free energy to describe the 
transition between a homogeneous and a periodically ordered system [l4l . [l5j . For a considerably earlier application 
to phase transitions in amphiphilic systems, see Ref. (l6| . Its central features are: 

• for €2,04 > 0, the term oc 0V 2 </> favors a periodically varying (f> and the term oc ^>V 4 ^> punishes a spatially 
varying <p 

• depending on the parameters, it may have as equilibrium states periodically ordered phases in one dimension 
(stripes), two dimensions (rods) and three dimensions (bec, fee, hep) 



• the characteristic wavenumber of the order parameter field is go = y^al (2C4), which follows from 

4>{v) (c 2 V 2 + C4V 4 ) 4>(v) = 4>{v) (- C4 g 4 + Ci {ql + V 2 ) 2 ) 0(r) 

It turns out that the phase diagram of the above model is equivalent to the formulation of a reduced model with the 
free energy according to [l7| 



^pfc = j d 3 xf PFC = J d 3 x (V(x) [-c + (1 + V) 2 ] (x) 



*(x) 



(18) 



which we call the (actual) phase field crystal (PFC) model. Indeed, we can define the dimcnsionless coordinate x, a 
free energy scale Eq and the reduced field ^ through the transformations 




E = , (19) 



9id 
t 



and the free energy in Eq. (fT?)) becomes 

(3T = E J d 3 x (B + Si*(x) + /pfc) • (20) 
with the value of e (see Eq. (|18jl for the definition of /pfc) given by 



The constants Bq and B\ are given by 



b 2 

aid - Po( c o - c A ql) + ) . (21) 



o 108gf d /?f + 36gf d b idPo (3p + 6 2 d (6a id gi d - fr 2 d 4- 6g id pgco) 

^0 = moj „8„4„2 > V-^j 



W8gt d q^c 
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o /""^ 4 27gf d p (3p + 6id(9a id 5id - 26? d + 9g id pgc ) , 
«i = yPod%9id 9? 2 8 4, • (23) 

The term Bo + in the free energy Eq. ([217)) does not influence the location of the phase boundaries since it is 
only linear in but it affects the values of the free energy density and the chemical potential at coexistence. Since 
we will determine these values explicitly lateron, we have given the expressions for Bo and B\ explicitly. Thus the 
phase diagram is determined by the variables occuring in /pfCj i-e. only by the parameter e and 'J, the average value 
of ^(x). The associated phase diagram has been calculated in Ref. [131 an d is depicted in Fig.[TJ 

In the next section, we turn to an exemplification of DFT, Taylor-expanded DFT and PFC for the hard sphere 
system. In particular, we will examine minimized solutions for crystals, values for crystal-liquid coexistence and 
crystal-liquid surface tensions as well as interfacial profiles for density modes. It will turn out that the apparent 
straightforwardness of the PFC derivation from Taylor-expanded DFT is misleading and cannot be upheld, if one 
likes to describe crystals with PFC. 
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FIG. 1. (color online) Phase diagram of the PFC free energy (|18[1 . Data are taken from Ref. [13]. The horizontal lines mark 
the values e = 0.5, 0.65 and 0.8 for which explicit results are discussed (see below). 



III. RESULTS FOR HARD SPHERE CRYSTALS AND CRYSTAL-LIQUID INTERFACE 



A. DFT: Fundamental Measure Theory 

For hard spheres, fundamental measure theory (FMT) allows the construction of very precise functional fl8rl2"ll ]. 
Essentially, FMT postulates an excess free energy with a local free energy density expressed in a set of weighted 
densities n a : 

^ cx [p] = j 'd 3 r$(n a (r)) . (24) 

The weighted densities are constructed as convolutions of the density with weight functions, n a (r) = p * w a (r). The 
weight functions reflect the geometric properties of the hard spheres. For one species, the weight functions include 
four scalar functions w°...w 3 , two vector functions w 1 , w 2 and a tensor function w 1 defined as 

3_/)/ D i„n „.,2_^ D i„n „..i _ w2 „.,o _ w2 



w a = 0(R-\y\), vf = 8{R-\r\), w L = — — , w 



4nR ' AttR 2 ' 



r 



<- = ^(^-H)- (25) 



Here, R = a/2 is the hard sphere radius. Using these weight functions, corresponding scalar weighted densities 
no . . .713, vector weighted densities ni,n.2 and one tensor weighted density n t arc defined. In constructing the free 
energy density $, arguments concerning the correlations in the bulk fluid and arguments for strongly inhomogeneous 
systems are used (for reviews see Refs. [2(J[2lj]). For the bulk, $ is required to reproduce exactly the second and third 
virial coefficent of the direct correlation function. Furthermore, consistency with a scaled particle argument and/or 
imposition of the Carnahan-Starling equation of state leads to the following form of the excess free energy density 

$({n[p(r)]}) = -n ln(l - n 3 ) + ^(n 3 ) " 



1 - n 3 

^2(^3) — ■ — ' 71 rr 1 — — — ■ (26) 



3 (—712 n 2 • n 2 + Ti2,int ,jjri2j + n 2 n t j jn t ji - rH,ijTk,jknt,ki) 
167r(l-n 3 ) 2 



Here, (^3) and ^2(^3) are functions of the local packing density 713 (r). With the choice 



(pi = 1 and ip2 = 1 



(27) 
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1.026 




0.94 
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0.65 


0.80 




3 

per cr 


1.039 


1.049 
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1.04 




1.041 [23J 1.039 [24] 


Pfier 3 


0.945 


0.944 


1.021 




0.94 




0.940 [23J 0.938 [24] 


(WO* 


4.96 


5.33 


7.23 




5.20 




4.96 f 




3.82 


3.99 


5.05 




3.93 




3.75 § 


/3/-icocx 


16.38 


17.44 


21.51 




17.16 




16.09 § 




2 • 10" 5 


0.10 


0.09 


-0.11 


-0.12 


-0.13 


3 • 10 -4 [25] 


/3o" 2 7[ioo] 


0.69 [26] 


0.89 


1.31 


0.140 


0.074 


0.046 


0.58 [24] 0.63 [26J 0.64 [27] 


/3o" 2 7[iio] 


0.67 


0.85 


1.21 


0.132 


0.070 


0.043 


0.56 0.61 


/3o" 2 7[iii] 


0.64 


0.78 


1.09 


0.105 


0.055 


0.034 


0.54 0.60 



t Free energy for pci-a" 3 = 1.041 [53] using an improved fit in the form of the Speedy equation of state from Ref. [s|. 
^ Free energy and chemical potential for pna 3 = 0.940 [23|] from the Carnahan-Starling equation of state. 

TABLE I. Results for coexistence properties and crystal-fluid surface tension using the different approaches considered in 
this work. The PFC results have been obtained using the fitting procedure described in Sec. MI CI for the reference density 
poo" 3 = 0.94 with the crystal and fluid coexistence densities as input (in italics). 



we obtain the Tarazona tensor functional [19( which is built upon the original Rosenfeld functional [l8| . The latter 
gives the fluid equation of state and pair structure of the Percus-Yevick approximation. Upon setting 

, , 2n 3 -n§ + 2(l-n 3 )ln(l-"3) (0R , 

<Pl - 1 H n \ Z °) 

3n 3 



^2 



2n 3 - 3n| + 2n\ + 2(1 - n 3 ) 2 ln(l - n 3 ) 



3n§ 



we obtain the tensor version of the recently introduced White Bear II (WBII) functional [22( . 

For bulk crystals, very accurate free energy results can already be obtained using a Gaussian approximation for the 
density, 

p(v) = E (1 " ^ac) (^) § exp (-«(r - vrf/a 2 ) , (29) 

lattice sites i 

and minimizing the total free energy with respect to the width parameter a and the vacancy concentration n vac at 
a fixed bulk density. At coexistence a ~ 80, in accordance with simulations and free energies per particle between 
Gaussian parametrized DFT and simulations agree on the level of 0.01 k B T @. However, the Tarazona functional ([27]) 
yields n vac — > while the WBII functional Eq. ([25)1 gives a finite equilibrium vacancy concentrations n vac = O(10 -5 ) 
which is about a factor 10 smaller than in the simulation results [25:]. This fine difference has important consequences: 
performing an unconstrained minimization (see Eq. (|10[) with vanishing external potential) only the WBII functional 
gives an absolute minimum for the free energy with a value for the chemical potential which is consistent with the 
derivative of the crystal free energy density with respect to the bulk density (see Ref. Q for further details). This 
implies that a free minimization for the crystal-fluid interface can only be performed using the WBII functional. 

The free minimization of the crystal-fluid interface is a non-trivial task, a brief description of the method (also 
applicable in the case of Taylor-expanded DFT) is given in App. Results for the surface tension with crystal 
faces in different orientations have been reported in Ref. (26[ and are also given in Table Q] There is agreement with 
simulations in the ordering 7[ioo] > 7[no] > 7[m] an d as far as the accuracy of the data permits, also in the values 
of the relative anisotropies (i.e. the values of (7[ioo] ~ 7[no])/7[ioo] an d (7[iool — 7[m])/7[ioo])- There is no clear 
consensus between different simulation methods on the absolute values of the 7's, the latest results are closer to the 
DFT values, however. 

Overall these results corroborate what is known from many applications of FMT on (dense) liquids [U, [28|, [29| : it 
is a quantitative theory and the accuracy also extends to the description of crystalline systems. Therefore we can 
consider FMT as a benchmark theory against which subsequent approximative approaches should be tested. 
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FIG. 2. (color online) Direct correlation function for hard spheres at a density per 3 = 0.9, somewhat below freezing. Simulation 
data are taken from Ref. [30l ]. 



B. Taylor-expanded DFT 

We consider the Taylor-expanded functional Eq. (fT2"j) of the FMT functional, Eq. ([27)1 (Tarazona functional) 
and Eq. (|28l (WBII functional). The nontrivial, second-order term in the functional involves the direct correlation 
function c^'{r\ p Q ). It is the second derivative of the excess free energy functional and is given in both cases by the 
polynomial form 



.,(2) 



(r; po) = (oi + a 2 r/a + a 3 (r/a) 3 )6(a - r) . (30) 



Using the packing fraction r/ — it a 3 / 6 po, the coefficients for of the Tarazona functional are those of the famous 
Percus-Yevick solution, 

(I + 277) 2 

ai — 



(i-v) 4 

1 + V/' 

(1 - v y 



(1 + ??/2) 2 

a 2 = 677— — r^— , (Tarazona) (31) 



V 

«3 = • 



The coefficients for of the WBII functional are given by 

1 + Arj + Arj 2 - 4r/ 3 + if 



ai 



03 = 



(1-77)4 

-2 + 2577 + 12?7 2 - IO77 3 + 2?7 4 2 ln(l - rf) 

3(1 - ?7) 4 + 3^7 

1 - 4r? + 2t7 2 - 3?7 3 + ?7 4 ln(l - 77) 



(White Bear II) (32) 



(I-??) 4 



Both forms approximate the simulation results for reasonably well, see Fig. ® with the WBII form matching 
more closely. As we will see, this does not imply better results for the Taylor-expanded functional in the WBII case. 

The Taylor-expanded functional contains the reference density po as an additional parameter. It would be desirable 
to choose it such that it also recovers the fluid density at coexistence. In this way it is guaranteed that at least the 
fluid properties are almost exact if the Taylor-expanded functional is fixed adequately. In order to determine po with 
moderate effort, we employ the approach already taken in Ref. (3lj . For a given bulk density p^, we parameterize 
the density profile in the Gaussian form Eq. (|29p and determine the three free parameters a (Gaussian width), n vac 
(vacancy concentration) and po by minimizing the functional for the grand potential difference per particle between 
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crystal and liquid at the chemical potential (J-(po) 



N 



1 f JS.f. _P(') 



d d r p(r) In p(r) + p 



J_ /" d 3 rdV C ( 2 )(r,r';p )Ap(r)Ap(r') 



(33) 



where Ap(r) = p(r) — p is defined as before. In practice, the space integrations can be restricted to one cubic fee 
unit cell. Finally, the bulk density pb is varied until AQ(pb) = 0, i.e. at this point pb and po (from the minimization) 
correspond to the coexisting crystal and fluid densities within the Gaussian approximation. We pick this po as 
the reference density for the unconstrained minimization, and re-determine the bulk crystal and fluid densities at 
coexistence through a full minimization and a subsequent Maxwell construction. These densities are only slightly 
shifted from the ones obtained within the Gaussian approximation (see Table [I]). 

The results are partially surprising (see Table HJ. For the Taylor-expanded Tarazona functional (i.e. employing 
the c^ 2 ) of Percus-Yevick) , the coexisting densities are still very close to the simulation and FMT (0.944 and 1.049). 
In the Gaussian approximation, similar values have been already obtained in 1985 by Jones and Mohanty (3lj . The 
crystal free energy is too big {f3F/N = 5.33 vs. 4.96 from simulation and FMT) and the width of the Gaussian peaks 
is much too narrow (a ~ 600 vs. 80 from simulation and FMT). For the WBII functional (i.e. the better functional, 
with a more precise c^) the coexistence densities are considerably off (1.021 and 1.123), consequently the crystal 
free energy is too big by 40% and the Gaussian width parameter a ~ 1000 stands for even narrower peaks. For both 
functionals, the vacancy concentrations are too large by 3 orders of magnitude. 

As before for FMT, we determine the surface tensions using a full minimization (see App. |A"|) . It is gratifying that 
the ordering 7noo] > 7[no] > 7[m] is upheld but the relative anisotropics are too large by approximately a factor of 2 
and the average value of the surface tension is too large (Taylor-expanded Tarazona functional: 0.84, Taylor-expanded 
White Bear II functonal: 1.2 vs. 0.67 from full FMT, all values in units of l/(/3er 2 )). 

In conclusion, it is apparent that through Taylor expansion of the precise FMT functionals, crystal free energies, 
surface tensions and vacancy concentrations are severely affected. Nevertheless, a qualitative descriptions is still 
achieved. The reason for the quantitative failure of the Taylor-expanded functionals is most likely due to the fact 
that the packing constraints or free energies for highly localized states are not respected ver y w ell. This is in contrast 
to the FMT functionals which have incorporated the correct description of localized states jl9l |20T ] . 



C. Phase-field crystal modelling 

1. Fixing parameters using bulk properties 
From the phase diagram of the PFC model (Fig. (fT])) with free energy density 

jp F C = *(X) [-€ + (1 + V) 2 ] Vf(x) + 

one infers that, for the description of stable fluid-crystal(fcc) coexistence, a parameter range of e ~ 0.5...1 is necessary. 
For lower e, fee is only metastable with respect to bec, and for e < there arc no ordered phases at all. However, 
following the derivation of PFC from the Taylor-expanded functional, one is led to the free energy in Eq. (fT7|) 
with coefficients aid, &id, Sid from the expansion of the ideal gas free energy and the Taylor-coefficients cq, c%, C4 
from the expansion of the Fourier transform of the direct correlation function c^ 2 '(k;p). Using the explicit form for 
c ( 2 )(r) = a,\ + a 2 (r/a) + a 3 (r/a) 3 (r/a < 1), we find 



(34) 



and inserting into the equation (f2"Tj) for e we obtain e = —0.6 ■ • • — 0.7 (po<r 3 = 0.94 . . . 1.04, Percus-Yevick direct 
correlation function). Thus, the naive gradient expansion of the Taylor-expanded functional produces a free energy 
which shows no sign of a liquid solid transition! The reason is essentially that the gradient expansion roughly 
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approximates the Fourier transform of the direct correlation function c( 2 \k; p) and consequently also the structure 
factor defined by 

S{k]P)= l-p^(k;p) - (35) 

The gradient expansion leads to a structure factor which clearly violates the empirical Verlet-Hanscn freezing criterion 
(height of first peak in S(k) > 2.8) whereas S(k) from the PY direct correlation function fulfills it. 

Such a failure of the naive gradient expansion has been noted and discussed before in a case study on the applicability 
of PFC for bcc metals The remedy proposed was to fit c , c 2 , c 4 to the first maximum of cf- 2 \k; p) (or S(k; p)) at 
around ka ~ 7 and also to fit the coefficients aid, &id, <?id m order to achieve a reasonable description of coexistence. 
The results from this procedure can be considered as partially successful: the correct description of the first peak 
of the structure factor needs a value for Co which is too small and hence causes deviations of the liquid isothermal 
compressibility ^(dp/dp)^ 1 = 1/(1 — pc^(0; p)) = S(0;p). Also, bulk free energies are not well captured [lol ]. 

Another serious problem is related to the identification of the order parameter <E>(r) = p(r)/po — 1 (see Eq. (jXTJ)) 
with the shifted dimcnsionlcss density. The order parameter ^ of the PFC model is related to a rescaling and shift 
of the order parameter </>. Numerical solutions of the PFC model for e > show that the order parameter solutions 
for bulk crystals can still be approximated by Gaussians (Eq. (|29[1 ). but they are much more spread out (width 
parameter a ~ 10, compared to 80 (for the case of FMT) and 500. . . 1000 (for the case of Taylor-expanded DFT). 
Consequently <&(r) has to be interpreted rather as a smeared-out reduced density. We will model this idea using a 
simple, normalized Gaussian smearing function with width a' leading to the following reinterpretation of $(r): 

$(r) = ^ - 1 , (36) 
Po 

with p(r) = J dr' 2 exp(-aV' 2 /CT 2 )p(r - r') . (37) 

Inserting this ansatz into the approximative free energy in Eq. (|17p and applying = —5J- c */(5p5p), we find the 
following Fourier transform for the direct correlation function: 

c& (k; p) = cxp (-^7) (-co + c 2 k 2 - c 4 k 4 ...). (38) 



Treating o/,co,C2,C4 as fitting parameters, we obtain very accurate fits of both, c^ 2 '(k;p) and S(k) in the "long 
wavelength" region ka < 10 for different choices of the reference density po, see Fig.|3]and Table [TT1 The good matching 
properties are a result of the fit with a' w 14 This value naturally accounts for the extended width of the bulk crystal 
solutions for $(r). Additionally we also have the qualitatively correct behavior for c^ 2 \r) for short distances r < a 
(see Fig. Furthermore, in line with the previous studies [13, H3 1 we treat the coefficients , , — > a,b,g as 
free parameters since the ideal gas free energy can not be applied to a smoothed density field p as defined above. A 
direct way to fix these parameters is by using the physical hard sphere coexistence densities p CI and pa as input. Since 
the PFC phase diagram is described by the reduced PFC free energy in Eq. ([15]). see Fig. [TJ the triple e, $f CI (e), ^a(e) 
fixes the three coefficients a,b,g using the third relations of Eq. (TT9")) and Eq. ([2~Tj) . Here, ^ cr (e) and ^fi(e) are the 
coexistence values for the average order parameter of the fee crystal and of the fluid, respectively. 
In summary, a reasonable procedure to fix the PFC parameters is as follows: 



fix a reference density po, fit S(k;po) = 1/(1 — poc^ 2 ' (k; po)) using Eq. (|38|) in the wave vector region including 
the first peak of the structure factor (ka < 10). 



• fix the PFC parameter e and require coexistence at the physical coexistence densities: this determines a, 6, g 
and consequently also the length scale go and the free energy scale Eq of the PFC model 

We have gathered the results of this procedure for some combinations of reference densities and e-parameters in 
Table HU Note that the fitted values of a, 6, g are one to two orders of magnitude larger than the ideal gas values 
«id = Po, = Po/2, ,9id = Po/3- This is a consequence of the fact that the PFC order parameter should be considered 
as a smeared-out density. 

In Table HI1 we also give the free energy per particle for the coexisting bulk liquid and crystal phases, their chemical 
potential and the vacancy concentration of the coexisting crystal, obtained from a minimization of the PFC free energy 
of Eq. (TT8"|) for bulk crystal states. . For the absolute values of the free energy and the coexisting chemical potential, 
one needs to determine the constant and linear terms in the PFC order parameter ^> (see Eq. (j2T))) ). We find that, for 
the reference density poa 3 = 0.94, the coexistence free energies and the chemical potential are well recovered in this 
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0.94 0.5 
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0.8 


14.21 61.93 2.567 0.0249 


23.45 69.75 33.13 0.279 
32.40 112.0 62.75 0.147 
42.02 160.5 101.0 0.091 


3.93 5.20 17.16 -0.11 
3.93 5.20 17.16 -0.12 
3.93 5.20 17.16 -0.13 


1.0 0.5 
0.65 
0.8 


13.93 80.36 3.201 0.03013 


21.04 98.19 54.57 0.342 
26.75 154.9 103.4 0.181 
32.42 218.6 166.3 0.112 


3.99 5.06 15.06 -0.06 

4.01 5.01 14.48 -0.07 

4.02 4.97 13.90 -0.08 


1.04 0.5 
0.65 
0.8 


13.99 96.56 3.745 0.03458 


16.58 122.9 76.15 0.397 
18.40 191.1 144.2 0.209 
19.53 266.5 232.1 0.130 


4.11 4.94 12.77 -0.03 
4.15 4.89 11.92 -0.04 
4.18 4.84 11.07 -0.05 



TABLE II. PFC fitting parameters and results for the bulk liquid and crystal phases at coexistence. 
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FIG. 4. Comparison between FMT and PFC for the order parameter in the bulk crystal at coexistence, for three different 
lattice directions i n the fee c ubic unit cell of side length a. The order parameter is calculated from the FMT density profile 
p(r) by *P(r) = ^/g/ (Po c 49o)(/^( r ) / P° ~ 1 — with q$ = C2/(2c4) and p(r) given by the convolution in Eq. ()37|l . Values 

for b, g, C2, ca are optained for the data set poo = 0.94 and e = 0.5 (see Table HI)) . and the smearing width for p is a 1 = 14.0 in 
(a) and a' = 7.0 in (b). 
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FIG. 5. (color online) Surface tension from the PFC model for different orientations and different reference densities po- For 
different po, the physical surface tensions 7 = -yqoEo differ since the free energy scale Eo and the inverse length scale squared 
Qo = C2/(2ci) differ (see Tab. [II}. The dimensionless PFC surface tension 7 is only a function of e (values given in text). 



case po = pa , the liquid free energy and the chemical potential are those of the PY theory and are hence reasonably 
accurate. The good approximation of the crystal-fluid free energy difference is in contrast to the findings in Ref. [1(| ■ 

However, the order parameter profiles $ from FMT and PFC in the crystal unit cell at coexistence agree only 
qualitatively, see Fig. [4] Here, the FMT solution for is calculated via the third expression of Eq. (|T9|) . but with 
the smeared density profile defined in Eq. (|37p . Using a smearing parameter a! = 14.0, consistent with the structure 
factor fit in Table ITl| yields order parameters that vary more strongly in FMT than in PFC (Fig. Eta)). Although in 
Fig. 2(a) the comparison is shown only for one choice of po and e, the differences between FMT and PFC are also 
found for other parameter combinations. Only by choosing a' = 7.0 (stronger smearing), the order parameter profiles 
agree almost quantitatively (see Fig. |UJb)). For such low values of a 1 , the fitting procedure gives too large values 
for the inverse PFC length scale go and too small values for the free energy scale Eq. Thus, the order parameter 
description in PFC is somewhat defective. 

The relative vacancy concentration n vac can be calculated via the number of particles A^ ce ii = /9 C r(amin/<7o) 3 in the 
fee cubic unit cell by n vac = 1 — N co n/4. Here, a m ; n is the cubic unit cell length in dimensionless PFC units which 
follows from minimizing the PFC free energy density at the value for the average order parameter ^ at coexistence. 
The resulting n vac is negative (see Table [H]) and of order 0.1 which implies a considerable concentration of interstitial 
particles. This is, of course, unphysical but it simply follows from fixing PFC coexistence to the correct physical 
densities. This observation reflects once more the difficulties in fixing parameters. Similarly to the case of the Taylor- 
expanded functional, it can not be expected that a "generic" free energy functional like the PFC functional can 
capture the correlations in the nearest-neighbor shell of a crystalline particle correctly, especially the condition of no 
overlap between particles. These correlations determine the precise value of n vac . 



2. Crystal-fluid surface tensions 



We have determined the equilibrium order parameter profile and the associated PFC free energy for the crystal-fluid 
interface as the long-time limit of the solution to the dynamic equation 

0*(x,t) _ 2 SF PFC 

(see Eqs. ([T]) and ©), with initial conditions given by a trial profile for a crystal slab in the simulation box filled 
otherwise with liquid. Any mobility coefficient, relating the PFC time t in the above equation to real time, is 
unimportant for the discussion of equilibrium properties. The order parameter profiles and the PFC free energy can 
be converted to density profiles and physical free energies using Tab. HU We have noted a certain sensitivity of the 
surface tensions and the order parameter profiles to the grid spacing and the precise extensions of the simulation box. 
These details are discussed in App. [B] 

We have calculated the dimensionless PFC surface tension 7 for the three different orientations [100], [110] and 
[111], each for the values of the PFC parameter e = 0.5, 0.65 and 0.8. In this sequence of orientations the results 
are: 0.0097, 0.0092, 0.0073 (e = 0.5), 0.0132, 0.0129, 0.0102 (e = 0.65) and 0.0165, 0.0163, 0.0129 (e = 0.8). The 
conversion to physical surface tensions of the hard sphere system is given by 7 = jq^Eo with the free energy scale Eq 
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and the inverse length scale squared q$ = 02/(204) given in Tab. HI1 This results in different 7 for different reference 
densities po which are depicted in Fig. [5j For the reference density po<J 3 = 0.94 (where PFC bulk crystal data are in 
good agreement with FMT and simulation) the surface tension values are given in Tab. Q] as well. The physical surface 
tensions are largest for highest reference density po and the lowest PFC parameter e and decrease for decreasing po 
and/or increasing e. Still, on average for the three different orientations the surface tensions are too low compared to 
FMT/simulation values between a factor of about 3 (po<r 3 = 1.04, e = 0.5) and about 15 (poa 3 — 0.94, e = 0.8). The 
ordering of surface tensions 7[xoo] > 7[no] > 7[m] i s correct for the PFC results but 7[m] is smaller than 7[ioo] by 
about 30% which differs considerably from the 5. ..8% as found in simulation and FMT. Likewise a strong qualitative 
difference between the order parameter profiles of the [111] interface compared to the [100], [110] interfaces is also 
found: the width of the [111] interface is considerably wider than of the [100], [110] interfaces (see below). This feature 
is not present in FMT. 

D. Density and order parameter modes at the crystal— fluid interface 

1. General theory and previous FMT results 

Consider a generic field ip(x,y,z) which describes the crystal-fluid interface with interface normal in z-direction. 
In DFT, this field is the density p(x, y, z) and in PFC, it is the order parameter field ^(x, y, z). We can parametrize 
the field ip in terms of a modified Fourier expansion 

i)(x, y,z)=Y^ exp(iKj • r) Pj (z) , (40) 

3 

where Kj denotes the reciprocal lattice vector (RLV), j and the z-dependent Fourier amplitude Pj{z) are modes of 
the field. One expects that upon crossing the interface from the crystal side, all Pj(z) relax to zero for nonzero Kj. 
Only for = 0, the value for the associated mode crosses from the average field ijj of the crystal to the average 
field of the fluid. It is convenient to group the K_, in shells with index m, where all Kj belonging to one shell, can 
be transformed into each other using the discrete symmetry group of the crystal under consideration. Thus in the 
bulk crystal, all Pj(z) = Pj associated with these Kj are equal. For the fee crystal, the reciprocal lattice is of bec 
symmetry We assume a to be the side-length of the cubic unit cell of fee and correspondingly 6 = 2ir/a the side 
length of the cubic unit cell of bec in reciprocal space. The reciprocal basis is given in Cartesian coordinates, where 
the axes span the cubic unit cell in reciprocal space, by Bx = 6(1, 1, —1), B 2 = 6(1, —1, 1) and B 3 = 6(— 1, 1, 1). An 
arbitrary RLV is a linear combination of the B^. The shells are characterized by a triple (to, n, k) of natural numbers 
and the Kj belonging to this shell have Cartesian components 6(±m, ±n, ±fc) and permutations thereof. Thus, if 
to, n, k are mutually distinct, there is a maximum of 48 RLV in one shell. The shells with lowest modulus are given by 
(1, 1, 1), (2, 0, 0) and (2, 2,0). A listing of the RLV triples up to shell 15 is given in Ref. ^ (Table I). At the interface, 
the degeneracy of the RLV in one shell is lifted, and we introduce an index n which distinguishes the possible values 
of the z-component of the RLV. Thus the decomposition becomes 

V, z) = / ]/,Pmn(z) exp (i(Kj) mn ■ r) . (41) 

mn j 

The sum over j is only for those RLV within shell to which have a common value of z-component, as expressed by the 
index n. In the literature, such a decomposition has been used to parameterize the full 3d density profile using only 
the leading mode in order to facilitate a simplified order parameter description of the crystal-fluid interface. In the 
context of PFC, the leading-mode picture has been advanced by Karma et al. [34| . If the z-component of (K.j) mn 
is zero, the mode will be purely real and if that z-component is nonzero, the mode will be in general complex and 
we denote by p^ in (z) its real part and by Pmn( z ) its imaginary part. The p mn (z) have the obvious interpretation of 
phase shifts of the associated field oscillations across the interface. 

The technique to perform the mode extraction from a full 3d solution il>(x, y, z) of a system with a solid-liquid 
interface is described in Ref. [35|- There, the density mode properties for the FMT solutions of the hard-sphere 
interfaces in different orientations have been discussed in detail. Some of these properties can be summarized as 
follows 

1. a separation of about one cubic unit cell length a ~ 1.6 a between the interface location as determined by the 
average density and the interface location as determined by the leading crystallinity mode (pi„(z)) 

2. a small density depiction zone just in front of the bulk crystal (dip in profile Poq(z)) 

3. strongly non-monotonic mode profiles also for next -to-leading modes, especially for P2n{z) 
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FIG. 6. Comparison between mode profiles from FMT and Taylor-expanded DFT (T-DFT, using the PY direct correlation 
function and poo~ 3 = 0.94). Note that the interface position for T-DFT is shifted by about 2a compared to the interface position 
for FMT to enhance readability. Real part of leading density modes for (a) the [100] interface and (b) the [111] interface. The 
density mode poo(z) has been rescaled and shifted, see tick labels at right t/-axis. Imaginary parts of leading density modes 
are displayed in (c) for the [100]-interface and in (d) for the [lll]-interface. 



4. kink position for higher modes p mn (z) shifts towards the bulk crystal for increasing m 



2. Mode profiles from T-DFT and PFC in comparison with FMT 

First, we compare the leading density modes of the [100] and [111] interfaces for the FMT solutions and the solutions 
from Taylor-expanded DFT (T-DFT), sec Fig.|6] The meaning of the different modes is seen best from the associated 
RLV. As a basis for the RLV, we choose vectors in direction of the cartesian axes with lengths 2n/a x: 2tt /a y and 2ir/a z , 
respectively. The quantities a x [y,z] are the side lengths of the minimal cuboid fee unit cells fitting the desired interface 
orientation in z-direction. For a graphical representation, we refer to Ref. [HI (Fig. 2). In Fig. [6j we illustrate the 
leading mode p\\{z) for the [100] orientation (corresponding to K = (1,1,1), the direction of close-packed planes) 
and the next-to-leading mode P22(z) (corresponding to K = (0,0,2), leading oscillation of lateral density average). 
For the [111] orientation, the leading mode splits into p\\ «->• K = (0,2,1) and pi2 K = (0,0,3). Both RLV 
correspond to directions of close packed planes, but the mode pti( z ) clearly differs from p^ 2 {z). Only the latter has 
a monotonic shape as expected for a "leading-order" interface profile, similar to P\i{z) of the [100] orientation. For 
the next-to-leading mode, we have P21 <-> K = (1, 1, 2). The FMT results for the real parts of these leading modes 
display already all properties 1-4 listed above. We further remark that the modes from T-DFT compare fairly well on 
a semi-quantitative level. The density depletion zone is missed and the phase shifts are more pronounced. For both, 
FMT and T-DFT, the mode expansion converges slowly as seen by the plateau values of the modes in the crystalline 
bulk. This is a consequence of the narrow density peaks in the bulk crystal. 

Next we compare the interface mode profiles for the order parameter vE'(r) from PFC on the one hand and from 
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FIG. 7. Comparison between order parameter mode profiles from FMT and PFC for the [100] interface in (a) and the [111] 
interface in (b) for the PFC parameter e = 0.5 (see Fig. [T]). The order parameter resulting from FMT is calculated from the 
3d FMT density profile by *(r) = y/g/(pfaq$) [p(r)/p - 1 - b/(3g)] (see Eqs. (O and {37]), ql = c 2 /(2c 4 )). We choose the 
reference density poa 3 = 0.94. The remaining parameters are given in Table [II] To enhance readability the "density" mode 
poo(z) is rescaled and shifted(see tick labels at right y-axis) and that the interface position for PFC is shifted by about 3a 
compared to the interface position for FMT. 

FMT on the other hand. Here, the FMT order parameter follows from the smeared density p(r) (Eq. (|37jl) which 
is rescaled and shifted according to the PFC transformation in Eq. (|19[) . For the PFC parameter, e = 0.5 and a 
reference density of poa 3 = 0.94, the results for the real part of the leading modes is shown in Fig. [7] (part (a) for the 
[100] interface and part (b) for the [111] interface). First one notes that the absolute magnitude of the PFC modes 
is smaller by a factor of about 3 compared to the FMT modes. This is a consequence of the different widths of the 
order parameter peaks in the crystal bulk: they are more narrow in FMT (see Fig. 0]) and consequently their Fourier 
amplitudes are larger. Secondly, except for the "density" depletion, the mode features identified in FMT are not 
present in PFC. This is attributed to the simple free energy of the PFC model. Our mode results further illustrate 
that the specificities of the layered hard sphere packing can not be captured by PFC. Another interesting observation 
is that the interface width w of the [100] interface is about 0.24 o for PFC and 0.86 a for FMT if the leading mode 
Pu(z) is fitted to the simple tanh-profilc 1 — tanh[(z — zo)/w]. For the [111] orientation we find widths of 0.47 a (PFC) 
and 1.05 a (FMT) from a fit to the leading mode pi2(z). In conclusion, we observe that the width of the interface is 
much smaller in PFC and it varies considerably with the orientation of the interface. 



3. Mode profiles from simulations in comparison with FMT 

We have carried out Molecular Dynamics simulations of the [100] interface in order to give a comparison to FMT 
data as well as to demonstrate the applicability of the mode expansion technique to simulation data. The simulations 
were carried out in the jVVT-ensemble at coexistence in cuboid boxes of cross-sectional area of 5 unit cells x 5 unit 
cells (L x x L y = 7.84cr x 7.84er) and a length of L z w 205<r with the crystal occupying about 60% of the box volume 
and placed in the middle of it. We have recorded the laterally averaged density profile 

P & y{z) = -j—j- / dx dyp(x,y,z) (42) 

with a resolution of 64 points per unit cell as a time average over different time intervals T av . 

From jOav(z) one can extract mode profiles p m n{z) for which the lateral components of the associated reciprocal 
lattice vector are zero: {K x ) mn — (K y ) mn — 0. In particular we focused on the average density mode poo(z) and the 
first two modes appearing for the lateral density average V2i{z) ((K) 2 2 = (0, 0, 2)) and y<o2{z) ((K) 6 2 = (0, 0, 4)). 

Due to the periodic boundaries of the simulation box, global centre of mass motion of the system does not cost any 
free energy, hence the system diffuses freely. When taking the time average, we have made no attempt to correct for 
this motion, as it was negligibly small on the time-scale of T av . Furthermore we have not corrected for the zero mode 
of the capillary waves at the interfaces. The zero mode corresponds to a shift Az of the average interfacial position, 
which is caused by fluctuations in the overall amount of crystalline material. For an infinite system, this zero mode 
would not incur a free energy penalty, either. For a finite system adsorption/desorption of crystalline layers results 
in a density change in the surrounding liquid reservoir. This is associated with a free energy cost, which we estimate 
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FIG. 8. Leading modes extracted from simulated, laterally averaged density profiles in comparison with FMT results, (a) 
average density mode poo(z), (b) real part of pin{z) ((K)a2 = (0,0,2)), (c) real part of p&2(z) ((K)22 = (0,0,4)) and (d) 
imaginary part of P2i(z). 



in quadratic approximation to be 

or j 

Ap = ^y {pa _ pfl) 2^ (pfl) (Az) 2 (43) 

where /ig (p) is the density derivative of the fluid chemical potential. For our system the broadening of the interface 
due to the zero mode contribution is rather small, ((Az) 2 ) < 1 a 2 . 

Besides the zero mode, there are also capillary waves with a finite wavelength. It is impossible to disentangle the 
contribution of capillary interfacial broadening from the width of a hypothetical "intrinsic" density profile which one 
would like to relate to the profile from density functional theory. However, by studying density averages for different 
time intervals, we will obtain some qualitative insight regarding the contributation that capillary waves make to the 
density modes. 

In Fig. |H]we show the mode profiles extracted from the simulation data for averaging times of 1, 5 and 50 self- 
diffusion times r as well as the FMT counterparts. (We give T av in units of the characteristic self-diffusion time, 
which it takes a particle in the coexisting liquid to diffuse over a distance of a.) The width of the average density 
(Fig. Ufa)) compares well with the FMT profile for T av = 1 and 5 but shows a significant broadening for T av = 50. 
This we largely ascribe to finite-wavelength capillary waves which are sampled better at longer times. The broadening 
effect on the mode width is less pronounced for the crystallinity modes (Fig.[5Jb)-(d)). A strong effect of the sampling 
time is visible on the plateau value in the crystalline part of the real parts of mode P22 (Fig. HJb)) and of mode ^62 
(Fig. Etc)). This reflects the broadening of the lattice site density peaks due to diffusion of the crystal as a whole. 
Apart from that the behavior of the modes in the interface region z/a — 7. . . 11 compares very well with the FMT 
results. In particular the good agreement for the imaginary part of P22 signifies that the wavelength shift of the 
density oscillations across the interface is captured correctly by FMT. 
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IV. SUMMARY AND CONCLUSION 



We have studied the crystal-liquid interface in the hard sphere system by means of theoretical description on three 
approximative levels, all of which are based on classical density functional theory (DFT): 

1. fundamental measure theory (FMT), the currently most accurate theoretical framework for hard spheres, 

2. second-order Taylor-expanded DFT which truncates density fluctuations with respect to a reference density 
beyond second order, 

3. a phase field crystal (PFC) model which formally can be viewed as local expansion of the second-order Taylor- 
expanded DFT in density fluctuations and in gradients thereof to fourth order. 

Coexistence properties, surface tensions and intcrfacial profiles for three interface orientations have been determined 
in all approximations by free minimization of the associated density functionals. FMT provides us with benchmark 
results against which the other approximations can be measured. For the [100] interface, we have confirmed that 
also the interfacial density mode profiles obtained from FMT show good agreement with corresponding results of 
Molecular Dynamics simulations. Thereby we have demonstrated the applicability of the mode expansion technique 
to simulation data. 

An important difference between FMT on the one hand, and Taylor-expanded DFT and PFC on the other hand 
is that the packing constraints for hard particles are incorporated very well only in the former. In consequence, the 
small values for the relative vacancy concentration for equilibrium crystals and the relatively small surface tension 
anisotropy are predicted correctly. Vacancy concentrations in Taylor-expanded DFT and PFC are off by orders of 
magnitude. The surface tension anisotropy is too large by a factor of two in Taylor-expanded DFT and much too 
large in PFC. In Taylor-expanded DFT, we find a surprising sensitivity of coexistence properties to details of the 
direct correlation function. This results in a strong sensitivity of the average surface tension with respect to the choice 
of the direct correlation function. 

We have discussed in detail the problem of parameter fixing in PFC. It turns out that the identification of the PFC 
order parameter with a rescaled and shifted smeared density is a suitable working recipe. The structure factor and 
coexistence free energies can be fitted very well. The order parameter distribution in the bulk crystal compares to 
FMT well only on a qualitative level, but with regard to the surface tension and the interfacial structure there are 
actually big discrepancies. Our major conclusion here is that the simple PFC variant considered here is too generic 
and needs to be specifically modified in order to incorporate the hard-sphere like interface structure of fee materials, 
(see Ref. [36| for an approach in this direction) . 

The sensitivity of results in Taylor-expanded DFT to the choice of the direct correlation function provides a hint 
that the specific functional form of the latter should be fitted to obtain proper coexistence. An attempt to match 
precise results for the direct correlation function (obtained by other means) is perhaps of little practical use. Wc further 
discuss that, at least for the hard sphere system, the status of Taylor-expanded DFT as the reference microscopic 
DFT for the PFC model is not justified. 

For future work on statics and dynamics of crystalline interfaces, grain boundaries and defects, it appears to be 
beneficial to pursue both, microscopically precise calculations using FMT and more coarse-grained investigations 
using the PFC model. Time and length scales accessible with FMT will be smaller than in PFC, but the PFC model 
can be better gauged in this way. 
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The free minimization of the crystal-fluid interface was achieved by an iterative solution of the Euler-Lagrange 
equation 



where the excess free energy J rcx is given by Eqs. (f2~i)) and (f2"6")l in the case of FMT and by Eq. (jT2"T) in the case of 
T-DFT, supplemented with the two choices for the direct correlation function (Eq. given by the coefficients 

in Eqs. (pITl) and (|3"2"1) . The density p(r) is discretized in a cuboid volume with edge lengths L x , L y and L z which 
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contains the fluid in the middle 

p{x,y,z) = p(x,y,L z - z). L x[l 



(z ~ L z /2) and the crystal phase at the boundaries (z ~ and z ~ L z ) such that 



are given by the edge lengths in :r[y]-direction, 



l x,[y] 



of the smallest cuboid unit 



cell of the crystal which has the desired orientation in z-direction. The crystal cuboid unit cells are shown in Fig. 2 
of Ref. [35| for the [100], [110], and [111] orientations. We chose L z = 32a z for the [100] and [110] orientations and 
L z = 16a z for the [111] orientation. The equidistant discretization was usually 64 points per unit cell length Ojuy, 
but was increased to 128 points per unit cell length a y \ z ] for the [111] case. All convolution integrals appearing in 
i5J rox /5p where calculated using 3d Fast Fourier transforms. 

As a first step for initialization, bulk crystal density profiles have been determined for the coexisting state. To do so, 
we determined the minimal free energy per particle of the bulk crystal in the cuboid unit cell at a fixed bulk density 
Ph by solving Eq. (|A1|) and also minimizing with respect to the unit cell length a (corresponding to a minimization 
with respect to the vacancy concentration [9(). The coexistence densities p cx , pa and the associated chemical potential 
M = Mcoex were determined using the Maxwell construction. The cuboid volume was filled with copies of the bulk 
crystal unit cell, defining a density profile /9 c (r). The initial interface profile po(r) was generated in the following way: 



p(z) = p a + (p CI - pa) p(z, z ,i,w{) , 
po{x, y, z) = p(z) + (p c (x, y, z) - p{z)) p(z, z , 2 , w 2 ) 
1 
2 



p(z,zo,w) 



(2 - tanh 


Z - Z 


— tanh 


L z — z — zq 


) 




w 




w 





(A2) 
(A3) 

(A4) 



The advantage of the above prescription is in the separation of the interfacial kink of the average density p from the 
kink of the density oscillations through different choices for zq,i and 2:0,2- This is useful to ensure a smooth start into 
the iterations in the case of FMT, for T-DFT it is not that important. 

Iteration was done using a combination of Picard steps with variabe mixing and DIIS steps (discrete inversion in 
iterative subspace) [37j- The Picard steps were performed according to 



p i+ i(x,y,z) = a{z) K[pi(x,y,z)] + (1 - a(z)) pi(x,y,z) 
a(z) = a m i n + (a ra ax - a m in)p'(z, z , w) , 



2 - tanh 2 


Z - Z 


- tanh 2 


L z — z — zq 




w 




w 



(A5) 
(A6) 

(A7) 

The mixing function a(z) ensures that there are substantial changes within one iteration step only in the interfacial 
region, since we chose a m in ~ 10 -5 and a max ~ 0.001. ..0.01. Choosing a(z) = const, is not practical since the constant 
would be limited to values below 10 -4 , otherwise the iteration fails due to instabilities in the bulk crystalline region. 
The DIIS steps were performed using between nous = 7... 10 previous profiles. 

A typical FMT run consisted of an initial Picard sequence with about 50 steps and a max = 0.01. Then, we 
alternated between Picard sequences of a minimum of 10 steps and one DIIS step (which needs another nous Picard 
initialization steps). We decreased a max after each switch from DIIS to Picard by a factor 1.5 until a max = 0.001 was 
reached. Also, we varied the maximum position of the mixing function a(z) by choosing Zq randomly in a certain 
interval (located in the interface region) with a width of about 2 a after each switch from DIIS to Picard. This was 
done to overcome being trapped at intermediate profiles where the surface tension 



1 



7 



dx 



dy 



dz (/rV(r)(lnp(r) - 1) + / cx [p(r)] - fip(r) + Pcocx ) 



(A8) 



hardly changes between iteration steps (/ cx is the excess free energy density and p C oox is the coexistence pressure). 
The condition for switching from Picard to DIIS was that after the minimum of 10 steps the convergence parameter 



1 



a x a y L z jq 



dx I dy dz {K\pi(r)\ - pi{v)Y 
>o Jo 



(A9) 



was decreasing between subsequent steps. If not, the Picard iterations were repeated with another 10 steps until that 
condition was met. Otherwise DIIS might take one away from the equilibrium solution easily. The DIIS step usually 
resulted in a very noticeable change in 7 and also in in the subsequent Picard steps. It was, however, not possible 
in general to perform a second DIIS step immediately after the first one since the density profile obtained after this 
DIIS step lead to singularities in the free energy (local packing fraction nz > 1). We stopped the run when ej < 10 -3 . 

We emphasize that only through the procedure outlined above we were able to determine equilibrium profiles for 
FMT. The standard method for solving DFT, simple Picard iterations with possibly variable, but spatially constant 
mixing a, simply fails. Also without DIIS we were not able to arrive at equilibrium profiles within a reasonable time. 

For T-DFT, the above procedure does not seem to be necessary but resulted in a very quick convergence. 
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Appendix B: Minimization of the PFC free energy for the crystal— fluid interface 



In PFC, we perform simulations with periodic boundary conditions in each direction, as we do in DFT. In the crystal 
phase, this implies that a stress will be acting on the crystal unless the dimensions L x ^ z j of the cuboid simulation 
box fit exactly multiples of the corresponding unit cell lengths of the equilibrium crystal. In order to avoid this stress, 
we use a simulation box which is minimizing the free energy of the crystal, i.e. we determine the minimizing length 
of the cubic unit cell a m - m := 2-K/q of the fee crystal (given in dimensionless PFC coordinates, x = go 1 ")- F° r a given 
average order parameter we apply Brent's method to find the box length which minimizes the free energy. 

A test for a single crystal cubic unit cell in [100] orientation (for e = 0.5 at coexistence, = —0.448336) with 
numbers of points per direction TV = 8, 16, 32 and 64 has shown that numerical box effects disappear for cubes of 



(16) 



0.539469, 



edge length 16 and larger. The results for the reciprocal lattice parameters q are = 0.539898, q 
q(32) _ 0.539476 and c^ 64 -* = 0.539468. It is interesting to compare these numbers to the corresponding numbers 
obtained by expanding the crystal order parameter in reciprocal lattice vectors (see Sec. IIII D II) and cutting the 
expansion at a maximum number n s h for the reciprocal lattice vector shells. We find for n s h = 4, 6, 8 and 10 the 
values = 0.53990, q^ = 0.53989, qW = 0.53956 and q<- 1Q '> = 0.53948. This demonstrates that n sh corresponds 
roughly to N/2 and that for precise numerical results the few-mode approximation is not quite sufficient. 

In order to avoid numerical artifacts, we determine the minimal reciprocal lattice parameters q separately for each 
orientation; we simulate one unit cell of the crystal with ./V = 32 for the [100] and the [110] interface. The crystal 
unit cell in [111] orientation is simulated in a box with discretization 32 x 64 x 64. The cuboid crystal unit cells for 
the different orientations are the same as used in the DFT calculations (see Fig. 2 of Ref. [35[). 

For the initialization of simulations of the crystal-liquid interface, half of the simulation box is filled with a one 
mode approximation of the crystal, in [100] orientation given by 

(Bl) 



<3>(x) = "Jcr + Acos(qx) cos(qy) cos(qz), 



and the liquid part in the other half has the constant value at coexistence \&a. The box length in the ^-direction 
(perpendicular to the interface) is 32 crystal unit cell lengths for the [100] and [110] orientation, resulting in a 
simulation box with a total number of points of 32 x 32 x 1024. For the [111] interface, we use 16 unit cell lengths 
resulting in a box with a total number of points of 32 x 64 x 1024. The crystal resides in one half of the box, so that 
one interfaces is in the middle of the box and the other near the periodic boundary. 

The PFC simulation evolves according to the dynamic equation ([39f until the system relaxes. As an indicator 
for the relaxation, we use the average deviation Sfj, of the local chemical potential /u(x) = SFp-pc/S^(x.) from the 
coexistence value fi cocx and stopped the computation when Sfi ~ 10 -4 . 

For the calculation of the dimensionless surface tension 7 we use the formula 
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d 3 a 



f - /. 



(B2) 



from [32|, Eq. (50)], where / is the PFC free energy density (/ cr for the crystal at coexistence and /a for the liquid at 
coexistence). ^ denotes the PFC order parameter with ^> CI the order parameter average in the coexisting crystal and 
$9 the corresponding average in the coexisting liquid, is the interface area (in dimensionless PFC units). Upon 
reordering Eq. (|B2j) we find: 



2 7 =i/^ 



/ + 



/cr^fl ~ /aj&cr g 



-h 



(/cr - /fi) 



(B3) 
(B4) 



define /, ^ as volume averages of the free energy density and order parameter 



1 
O 

V 

n 



where V = J d 3 x = L x ■ L y ■ L z 



I 



^cr - * fi 
/crjfl ~ /ajger _ Jc 

L r , and we obtain 



/cr ~ /fl 
*cr - *fl 
ffi 



(B5) 
(B6) 
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TV 


7[ioo] 


8 


0.00913 


16 


0.01041 


32 


0.01041 


64 0.01052 



TABLE III. Surface tension 7 for e = 0.53 in [100]-direction. 



Note that the factor \ is needed due to the presence of two interfaces in the simulation. 



To calculate the surface tension with Eq. (|B7j) we calculate /, "J in the whole domain. ^>a, ^> cl -, f CI and fa are 
calculated by convoluting / and ^ with normalized Gaussians of sufficient width such that the resulting profile is 
locally constant. This is equivalent to peak to peak averaging of the / and ^ profiles on the crystal side. 



For e = 0.53, results for the surface tension have been reported previously in Ref. {3a ]. We checked the convergence 
of the surface tension for different number of points N per unit cell len gth . For the [100] orientation, the results are 
given in Tab. lIIII and should be compared with 7[ioo] = 0.0113 from Ref. |38j . For the [11 Indirection Ref. [IU provides 
7[m] = 0-0082 whereas our result is 7[m] = 0.0079 (using N = 32). It is not clear which precise discretization was 
used in Ref. [38| but we can conclude that typical discretizations of about 10. .15 points per unit cell used in the PFC 
community leave a residual error of about 5 per cent in the value of the surface tension. 
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